# Turn on the gmri font for plots
showtext::showtext_auto()
This report was created to track the sea surface temperature regimes for marine regions of interest to the Gulf of Maine Research Institute. The default region being a central snapshot of the Gulf of Maine.
Satellite sea surface temperature data used was obtained from the National Center for Environmental Information (NCEI). With all maps and figures displaying NOAA’s Optimum Interpolation Sea Surface Temperature Data.
The spatial extent for Gulf Of Maine is displayed below. This bounding box is the same bounding box coordinates used to clip the OISST data when constructing the time series data from the array.
# File paths for various extents based on params$region
region_paths <- get_timeseries_paths(region_group = "gmri_sst_focal_areas", mac_os = "mojave")
# Load the bounding box for Andy's GOM to show they align
poly_path <- region_paths[[params$region]][["shape_path"]]
region_extent <- st_read(poly_path, quiet = TRUE)
# Pull extents for the region to set crop extent
crop_x <- st_bbox(region_extent)[c(1,3)]
crop_y <- st_bbox(region_extent)[c(2,4)]
# Expand the area out to see the larger patterns
crop_x <- crop_x + c(-2, 2)
crop_y <- crop_y + c(-0.75, 0.75)
# Zoom out for cpr extent
if(tolower(params$region) == "cpr gulf of maine"){
crop_x <- c(-70.875, -65.375)
crop_y <- c(40.375, 45.125)}
# Add the bottom contours:
bathy <- raster("~/Documents/Repositories/Points_and_contours/NEShelf_Etopo1_bathy.tiff")
# Contours for geom_contour()
bathy_df <- as.data.frame(raster::coordinates(bathy))
bathy_df$depth <- raster::extract(bathy, bathy_df)
bathy_df$depth <- bathy_df$depth * -1
contours_make <- c(50, 100, 250)
# Full plot
ggplot() +
geom_sf(data = new_england, fill = "gray90", size = .25) +
geom_sf(data = canada, fill = "gray90", size = .25) +
geom_sf(data = greenland, fill = "gray90", size = .25) +
geom_sf(data = region_extent,
color = gmri_cols("gmri blue"),
fill = gmri_cols("gmri blue"), alpha = 0.2, linetype = 2, size = 0.5) +
geom_contour(data = bathy_df, aes(x, y, z = depth),
breaks = contours_make,
color = "gray80") +
map_theme +
coord_sf(xlim = crop_x,
ylim = crop_y, expand = F)
Area-specific time series are the most basic building block for relaying temporal trends. For any desired area (represented by a spatial polygon) we can generate a time series table of the mean sea surface temperature within that area for each day. Additionally, we can compare how observed temperatures correspond with the expected conditions based on a climatology using a specified reference period.
# Use {gmRi} instead to load timeseries to tye up loose ends
timeseries_path <- region_paths[[params$region]][["timeseries_path"]]
region_timeseries <- read_csv(timeseries_path, col_types = cols(), guess_max = 1e6)
# format dates
region_timeseries <- region_timeseries %>%
mutate(time = as.Date(time))
# Display Table of last 6 entries
tail(region_timeseries) %>%
mutate(across(where(is.numeric), round, 2)) %>%
select(
Date = time,
`Sea Surface Temperature` = sst,
#`Area-Weighted SST` = area_wtd_sst,
`Day of Year` = modified_ordinal_day,
`Climate Avg.` = sst_clim,
#`Area-weighted Climate` = area_wtd_clim,
`Temperature Anomaly` = sst_anom#,
#`Area-Weighted Anomaly` = area_wtd_anom
) %>% gt() %>%
tab_header(
title = md(paste0("**", tidy_name, " - Regional Sea Surface Temperature", "**")),
subtitle = paste("Temperature Unit: Celsius")) %>%
tab_source_note(
source_note = md("*Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.*") ) %>%
tab_source_note(md("*Climatology Reference Period: 1982-2011.*"))
| Gulf Of Maine - Regional Sea Surface Temperature | ||||
|---|---|---|---|---|
| Temperature Unit: Celsius | ||||
| Date | Sea Surface Temperature | Day of Year | Climate Avg. | Temperature Anomaly |
| 2022-01-26 | 6.70 | 26 | 5.46 | 1.24 |
| 2022-01-27 | 6.19 | 27 | 5.35 | 0.84 |
| 2022-01-28 | 6.51 | 28 | 5.32 | 1.19 |
| 2022-01-29 | 6.62 | 29 | 5.25 | 1.36 |
| 2022-01-30 | 6.30 | 30 | 5.21 | 1.09 |
| 2022-01-31 | 6.28 | 31 | 5.18 | 1.10 |
| Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data. | ||||
| Climatology Reference Period: 1982-2011. | ||||
# march 1st sst
mar1 <- region_timeseries %>%
filter(modified_ordinal_day == 61) %>%
distinct(sst_clim) %>%
pull(sst_clim)
Each of our Climatologies are currently set up to calculate daily averages on a modified year day, such that every March 1st and all days after fall on the same day, regardless of whether it is a leap year or not.
This preserves comparisons across calendar dates such-as: “The average temperature on march 1st is 4.2378197` for the reference period 1982 to 2011”
In these tables Sea Surface Temperature is the mean temperature observed for that date averaged across all cells within the area. Climate Avg. & Climate SD are the climate means and standard deviations for a 1982-2011 climatology. Temperature Anomaly is the daily observed sea surface temperature minus the climate mean.
Regional warming trends below were calculated using all the available data for complete years beginning with 1982 through the end of 2021. The overlaid trend lines then track how warming has increased with time. A dotted line has been included to show how the global average temperature has changed during the same period.
# Summarize by year to return mean annual anomalies and variance
annual_summary <- region_timeseries %>%
mutate(year = year(time)) %>%
filter(year %in% c(1982:2021)) %>%
group_by(year) %>%
summarise(sst = mean(sst, na.rm = T),
sst_anom = mean(sst_anom, na.rm = T),
area_wtd_sst = mean(area_wtd_sst),
area_wtd_anom = mean(area_wtd_anom),
.groups = "drop") %>%
mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")))
# # Global Temperature Anomaly Rates
global_anoms <- read_csv(
paste0(oisst_path, "global_timeseries/global_anoms_1982to2011.csv"),
guess_max = 1e6,
col_types = cols()) %>%
mutate(year = year(time))
# summarize by year again
global_summary <- global_anoms %>%
group_by(year) %>%
summarise(sst = mean(sst, na.rm = T),
sst_anom = mean(sst_anom),
area_wtd_sst = mean(area_wtd_sst),
area_wtd_anom = mean(area_wtd_anom),
.groups = "drop") %>%
mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")))
# Build Regression Equation Labels
# 1. All years
lm_all <- lm(area_wtd_sst ~ year,
data = filter(annual_summary, year %in% c(1982:2021))) %>%
coef() %>%
round(3)
# 2. Last 15 years
lm_15 <- lm(area_wtd_sst ~ year,
data = filter(annual_summary, year %in% c(2007:2021))) %>%
coef() %>%
round(3)
# 3. Global - Area-weighted
lm_global <- lm(area_wtd_sst ~ year,
data = filter(global_summary, year %in% c(1982:2021))) %>%
coef() %>%
round(3)
# Convert yearly rate to decadal
decade_all <- lm_all['year'] * 10
decade_15 <- lm_15['year'] * 10
decade_global <- lm_global["year"] * 10
# Equation to paste in
eq_all <- paste0(decade_all, "\u00b0", "C / Decade")
eq_15 <- paste0(decade_15, "\u00b0", "C / Decade")
eq_global <- paste0(decade_global, "\u00b0", "C / Decade")
# Generate a smoothed temperature line using splines
yearly_temp_smooth <- as.data.frame(spline(annual_summary$yr_as_dtime, annual_summary$area_wtd_anom)) %>%
mutate(x = as.Date(x, origin = "1970-01-01"))
#### Annual Trend Plot ####
ggplot(data = annual_summary, aes(yr_as_dtime, area_wtd_anom)) +
# Add daily data
geom_line(data = region_timeseries,
aes(time, area_wtd_anom, color = "Daily Temperatures")) +
# Overlay yearly means
#geom_line(data = yearly_temp_smooth, aes(x, y, color = "Average Yearly Temperature"), alpha = 0.7, linetype = 2) +
geom_line(color = "gray10", size = 1) +
geom_point(color = "gray10", alpha = 0.7, size = 0.75) +
# Add regression lines
geom_textsmooth(data = filter(global_summary, year <= 2021),
method = "lm", text_smoothing = 30,
label = "Global Trend",
color = gmri_cols("green"),
linewidth = 1,
formula = y ~ x, se = F,
linetype = 3, hjust = 0.925) +
geom_smooth(data = filter(annual_summary, year <= 2021),
method = "lm",
aes(color = "1982-2021 Regional Trend"), #label = "40-Year Trend",
formula = y ~ x, se = F,
linetype = 2) +
geom_smooth(data = filter(annual_summary, year %in% c(2007:2021)),
method = "lm",
aes(color = "2007-2021 Regional Trend"),
formula = y ~ x, se = F,
linetype = 2) +
# Manually add equations so they show yearly not daily coeff
geom_text(data = data.frame(),
aes(label = eq_all, x = min(region_timeseries$time), y = Inf),
hjust = 0, vjust = 2, color = gmri_cols("gmri blue")) +
geom_text(data = data.frame(),
aes(label = eq_15, x = min(region_timeseries$time), y = Inf),
hjust = 0, vjust = 3.5, color = gmri_cols("orange")) +
geom_text(data = data.frame(),
aes(label = eq_global, x = min(region_timeseries$time), y = Inf),
hjust = 0, vjust = 5, color = gmri_cols("green")) +
# Colors
scale_color_manual(values = c(
"1982-2021 Regional Trend" = as.character(gmri_cols("gmri blue")),
"2007-2021 Regional Trend" = as.character(gmri_cols("orange")),
#"1982-2021 Global Trend" = as.character(gmri_cols("green")),
"Average Yearly Temperature" = "gray10",
"Daily Temperatures" = "gray90")) +
# Axes
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
# theme
theme(legend.title = element_blank(),
legend.position = c(0.825, 0.15),
legend.background = element_rect(fill = "transparent"),
legend.key = element_rect(fill = "transparent", color = "transparent"),
panel.grid = element_blank()) +
# labels + theme
labs(x = "",
y = "Sea Surface Temperature Anomaly",
caption = paste0("Anomalies calculated using 1982-2011 reference period."))
Without Daily Temperatures
#### Annual Trend Plot ####
ggplot(data = annual_summary, aes(yr_as_dtime, area_wtd_anom)) +
# Overlay yearly means and lines connecting them
geom_col(aes(fill = "Yearly Anomaly"), alpha = 0.7, size = 0.75) +
scale_fill_manual(values = c("Yearly Anomaly" = "gray90")) +
# Add regression lines
geom_textsmooth(data = filter(global_summary, year <= 2021),
method = "lm", text_smoothing = 30,
size = 4,
label = "Global Trend",
color = gmri_cols("green"),
linewidth = 1,
formula = y ~ x, se = F,
linetype = 1, hjust = 0.925) +
geom_textsmooth(data = filter(annual_summary, year <= 2021),
method = "lm",
# aes(color = "1982-2021 Regional Trend"), #label = "40-Year Trend",
size = 4, linewidth = 1,
label = "1982-2021 Regional Trend",
formula = y ~ x, se = F,
linetype = 1, linewidth = 1, color = gmri_cols("gmri blue"),
hjust = 0.50, vjust = -1.2
) +
geom_textsmooth(data = filter(annual_summary, year %in% c(2007:2021)),
method = "lm",
# aes(color = "1982-2021 Regional Trend"), #label = "40-Year Trend",
size = 4, linewidth = 1,
label = "2007-2021 Regional Trend",
formula = y ~ x, se = F,
linetype = 1, linewidth = 1, color = gmri_cols("orange"),
#boxlinetype = "dotted", boxlinewidth = 0.5,
hjust = 0.95, vjust = -1.2
) +
# Manually add equations so they show yearly not daily coeff
geom_text(data = data.frame(),
aes(label = eq_all, x = min(region_timeseries$time), y = Inf),
hjust = 0, vjust = 2, color = gmri_cols("gmri blue")) +
geom_text(data = data.frame(),
aes(label = eq_15, x = min(region_timeseries$time), y = Inf),
hjust = 0, vjust = 3.5, color = gmri_cols("orange")) +
geom_text(data = data.frame(),
aes(label = eq_global, x = min(region_timeseries$time), y = Inf),
hjust = 0, vjust = 5, color = gmri_cols("green")) +
# Axes
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
# labels + theme
labs(title = "Gulf of Maine Warming Faster then Global Average",
x = "",
y = "Sea Surface Temperature Anomaly",
caption = paste0("Anomalies calculated using 1982-2011 reference period.")) +
# theme
theme(legend.title = element_blank(),
legend.position = c(0.85, 0.1),
legend.background = element_rect(fill = "transparent"),
legend.key = element_rect(fill = "transparent", color = "transparent"),
panel.grid = element_blank())
# Lazy seasons
# season = factor(quarter(time, fiscal_start = 1)),
# season = fct_recode(season,
# c("Jan 1 - March 31" = "1"),
# c("Apr 1 - Jun 30" = "2"),
# c("Jul 1 - Sep 30" = "3"),
# c("Oct 1 - Dec 31" = "4")),
# Doing seasons by meteorological Definitions
quarter_summary <- region_timeseries %>%
mutate(
yr = year(time),
month_num = month(time),
month = month(time, label = T, abbr = T),
season = metR::season(month_num, lang = "en"),
season_eng = case_when(
season == "SON" ~ "Fall",
season == "DJF" ~ "Winter",
season == "MAM" ~ "Spring",
season == "JJA" ~ "Summer"),
season_eng = factor(season_eng, levels = c("Winter", "Spring", "Summer", "Fall")),
#Set up correct year for winters, they carry across into next year
season_yr = ifelse((season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr),
year = season_yr) %>%
filter(year >= 1982) %>%
group_by(year, season_eng) %>%
summarise(sst = mean(sst, na.rm = T),
sst_anom = mean(sst_anom, na.rm = T),
area_wtd_sst = mean(area_wtd_sst, na.rm = T),
area_wtd_anom = mean(area_wtd_anom, na.rm = T),
.groups = "drop")
# Plot
quarter_summary %>%
ggplot(aes(year, area_wtd_anom)) +
geom_line(group = 1, color = "gray60", linetype = 3) +
geom_point(size = 0.75) +
geom_smooth(method = "lm",
aes(color = "Regional Trend"),
formula = y ~ x, se = F, linetype = 1) +
stat_poly_eq(formula = y ~ x,
color = gmri_cols("gmri blue"),
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = T) +
scale_color_manual(values = c("Regional Trend" = as.character(gmri_cols("orange")))) +
labs(x = "",
y = "Sea Surface Temperature Anomaly",
caption = "Regression coefficients reflect annual change in sea surface temperature.") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
theme(legend.title = element_blank(),
legend.position = c(0.925, 0.05),
legend.background = element_rect(fill = "transparent"),
legend.key = element_rect(fill = "transparent", color = "transparent")) +
facet_wrap(~season_eng, ncol = 4)
dat_region <- annual_summary %>%
filter(year %in% c(1982, 2021))
dat_global <- global_summary %>%
filter(year %in% c(1982, 2021))
dat_list <- list(dat_region, dat_global) %>% setNames(c(tidy_name, "Global Oceans"))
dat_combined <- bind_rows(dat_list, .id = "Area") %>%
select(Area, area_wtd_sst, year) %>%
pivot_wider(names_from = year, values_from = area_wtd_sst)
ggplot(dat_combined, aes(x = `1982`, xend = `2021`, y = fct_rev(Area))) +
geom_dumbbell(colour = "lightblue",
colour_xend = gmri_cols("gmri blue"),
size = 3,
dot_guide = TRUE,
dot_guide_size = 0.5) +
labs(x = "Sea Surface Temperature",
subtitle = "Change in Sea Surface Temeprature - 1982-2021",
y = "") +
scale_x_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "temperature"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C"))
For the figures below heatwave events were determined using the methods of Hobday et al. 2016 and implemented using the R package {heatwaveR}.
A marine heatwave is defined as a situation when seawater temperatures exceeds a seasonally-varying threshold (usually the 90th percentile) for at least 5 consecutive days. Successive heatwaves with gaps of 2 days or less are considered part of the same event. The heatwave threshold used below was 90%. The heatwave history for Gulf Of Maine is displayed below:
# Use function to process heatwave data for plotting
region_heatwaves <- pull_heatwave_events(region_timeseries, threshold = 90)
For anything we wish to host on the web there is an option to display tables and graphs that are interactive. Interactivity allows users to pan, zoom, and highlight discrete observations.
# # Option 1: Grab last 365 days
#
# # Grab data from the most recent year through present day to plot
# last_year <- Sys.Date() - 365 #1 year from current date
# last_yr_heatwaves <- region_heatwaves %>% filter(time >= last_year)
# last_yr <- year(last_yr)
# # Option 2: Grab last Full Year:
# last_year <- Sys.Date() - 365 #1 year from current date
#
# # wind it back to first day of the year
# last_year <- last_year - yday(last_year) + 1
#
# # Filter out:
# last_yr_heatwaves <- region_heatwaves %>% filter(year(time) == year(last_year))
# last_yr <- year(last_yr)
# Option 3: Last complete year
last_year <- 2021
last_yr_heatwaves <- region_heatwaves %>% filter(year(time) == last_year)
# Get number of heatwave events and total heatwave days for last year
# How many heatwave events:
# How many heatwave events
num_events <- max(last_yr_heatwaves$mhw_event_no, na.rm = T) - min(last_yr_heatwaves$mhw_event_no, na.rm = T)
# Data for only the current year
this_yr_hw <- region_heatwaves %>% filter(year(time) == year(Sys.Date()))
# Number of heatwave events
num_hw_days <- sum(this_yr_hw$mhw_event, na.rm = T)
# Plot the interactive timeseries
last_yr_heatwaves %>%
plotly_mhw_plots()
# Set colors by name
color_vals <- c(
"Sea Surface Temperature" = "royalblue",
"Heatwave Event" = "darkred",
"MHW Threshold" = "coral3",
"Daily Climatology" = "gray30")
# Set the label with degree symbol
ylab <- "Sea Surface Temperature"
# Plot the last 365 days
hw_temp_p <- ggplot(last_yr_heatwaves, aes(x = time)) +
geom_segment(aes(x = time, xend = time, y = seas, yend = sst),
color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe),
color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
geom_line(aes(y = hwe, color = "Heatwave Event")) +
geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
scale_color_manual(values = color_vals) +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "temperature"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
theme(legend.title = element_blank(),
legend.position = "none") +
labs(x = "",
y = ylab,
title = str_c(tidy_name, " - ", last_year))
# Show Figure
# hw_temp_p
# Plot the last 365 days - anomaly scale
linetype_key <- c(
"Sea Surface Temperature Anomaly" = 1,
"Heatwave Event" = 1,
"MHW Threshold" = 3,
"Daily Climatology" = 2)
# Same Plot as Anomalies
hw_anom_p <- last_yr_heatwaves %>%
mutate(
sst = sst,
seas = seas,
sst_anom = sst_anom,
mhw_thresh = mhw_thresh,
anom_thresh = mhw_thresh - seas,
anom_hwe = hwe - seas) %>%
ggplot(aes(x = time)) +
geom_segment(aes(x = time, xend = time,
y = 0, yend = sst_anom),
color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time,
y = anom_thresh, yend = anom_hwe),
color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst_anom, color = "Sea Surface Temperature Anomaly")) +
geom_line(aes(y = anom_hwe, color = "Heatwave Event")) +
geom_line(aes(y = anom_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
geom_line(aes(y = 0, color = "Daily Climatology"), lty = 2, size = .5) +
scale_color_manual(values = color_vals) +
scale_linetype_manual(values = linetype_key, guide = "none") +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
guides(color = guide_legend(override.aes = list(linetype = linetype_key), nrow = 2)) +
theme(legend.title = element_blank(),
legend.position = "top") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
labs(x = "",
y = "Temperature Anomaly",
caption = paste0("Climate reference period : 1982-2011"))
# Show figure
# hw_anom_p
hw_temp_p / hw_anom_p
# Prep the legend title
guide_lab <- "Sea Surface Temperature Anomaly \u00b0C"
# Set new axis dimensions, y = year, x = day within year
# use a flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- region_heatwaves %>%
mutate(year = year(time),
yday = yday(time),
flat_date = as.Date(yday-1, origin = base_date))
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(grid_data$sst_anom)) *c(-1, 1)
# Assemble heatmap plot
heatwave_heatmap <- ggplot(grid_data, aes(x = flat_date, y = year)) +
# background box fill for missing dates
geom_rect(xmin = base_date, xmax = base_date + 365,
ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5,
fill = "gray75", color = "transparent") +
# tile for sst colors
geom_tile(aes(fill = sst_anom)) +
# points for heatwave events
geom_point(data = filter(grid_data, mhw_event == TRUE),
aes(x = flat_date, y = year), size = .25) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
labs(x = "",
y = "",
"\nClimate reference period : 1982-2011") +
#scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
scale_fill_distiller(palette = "RdYlBu", na.value = "transparent",
limit = limit) +
#5 inches is default rmarkdown height for barheight
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "right",
title.hjust = 0.5,
barheight = unit(4.8, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
theme(legend.title = element_text(angle = 90),
axis.line.x = element_blank(),
axis.line.y = element_blank())
# Assemble pieces
heatwave_heatmap
#### Annual Heatwave Summary Details
# number of heatwaves
# average heatwave duration
# remove NA as a distinct heatwave number
wave_summary <- grid_data %>%
group_by(year(time), mhw_event_no) %>%
summarise(total_days = sum(mhw_event, na.rm = T),
avg_anom = mean(sst_anom, na.rm = T),
peak_anom = max(sst_anom, na.rm = T),
.groups = "drop") %>%
drop_na() %>%
group_by(`year(time)`) %>%
summarise(num_waves = n_distinct(mhw_event_no),
avg_length = mean(total_days, na.rm = T),
avg_intensity = mean(avg_anom, na.rm = T),
peak_intensity = max(peak_anom, na.rm = T),
.groups = "drop") %>%
rename(year = `year(time)`)
#### Plotting
# number of heatwaves
hw_counts <- ggplot(wave_summary, aes(y = year, x = num_waves)) +
geom_segment(aes(yend = year, xend = 0),
color = gmri_cols("gmri blue")) +
geom_point(color = gmri_cols("gmri blue")) +
labs(x = "Event Count", y = "")
# average duration
hw_lengths <- ggplot(wave_summary, aes(y = year, x = avg_length)) +
geom_segment(aes(yend = year, xend = 0),
color = gmri_cols("orange")) +
geom_point(color = gmri_cols("orange")) +
labs(x = "Event Duration", y = "")
# avg temp
hw_temps <- ggplot(wave_summary, aes(y = year, x = avg_intensity)) +
geom_segment(aes(yend = year, xend = 0),
color = gmri_cols("green")) +
geom_point(color = gmri_cols("green")) +
labs(x = "Avg Anomaly", y = "")
# peak temp
hw_peaks <- ggplot(wave_summary, aes(y = year, x = peak_intensity)) +
geom_segment(aes(yend = year, xend = 0),
color = gmri_cols("teal")) +
geom_point(color = gmri_cols("teal")) +
labs(x = "Peak Anomaly", y = "")
(hw_counts | hw_lengths | hw_temps | hw_peaks) + plot_annotation(title = "Heatwaves Frequent and Powerful in Recent Years")
The 2021 global sea surface temperature anomalies have been loaded and displayed below to visualize how different areas of the ocean experience swings in temperature.
# Access information to netcdf on box
nc_year <- "2021"
anom_path <- str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", nc_year, ".nc")
# Load 2021 as stack
anoms_2021 <- stack(anom_path)
# Make Annual Average
ann_anom_ras <- calc(anoms_2021, fun = mean, na.rm = T)
# Color limit for palette
#temp_limits <- c(-5, 5)
temp_limits <- max(abs(values(ann_anom_ras)), na.rm = T) * c(-1,1) # Dynamic Limits
sst_lab <- expression("Sea Surface Temperature Anomaly"~~degree~C)
# Build Map
ggplot() +
geom_stars(data = st_as_stars(rotate(ann_anom_ras))) +
geom_sf(data = world, fill = "gray30", color = "white", size = .25) +
scale_fill_distiller(palette = "RdBu", na.value = "transparent",
limit = temp_limits, oob = scales::squish
) +
guides("fill" = guide_colorbar(title = sst_lab,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(3, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
coord_sf(expand = F) +
map_theme +
labs(title = str_c("Global Temperature Anomalies - ", nc_year))
#### Making Averages for Seasons
# Get previous year winter
last_nc_yr <- as.numeric(nc_year) - 1
last_yr_anom_path <- str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", last_nc_yr, ".nc")
# Load 2020 as stack to get the full winter
last_yr_anoms <- stack(last_yr_anom_path)
# Join to 2021
anoms_double <- stack(list(last_yr_anoms, anoms_2021))
# Drop december 2021 so its not influencing the Winter of 2020-2021
dec_2021 <- "X2021.12"
not_dec21 <- which(str_sub(names(anoms_double), 1, 8) != dec_2021)
anoms_nodec <- anoms_double[[not_dec21]]
# Set up list to use map()
month_key <- list("Winter" = c("12", "01", "02"),
"Spring" = c("03", "04", "05"),
"Summer" = c("06", "07", "08"),
"Fall" = c("09", "10", "11"))
# Get mean anoms across seasons
season_stacks <- map(month_key, function(mon){
# Get names from the stack
stack_names <- names(anoms_nodec)
stack_months <- str_sub(stack_names, 7,8)
# layers with correct months:
in_season <- which(stack_months %in% mon)
# Get mean across those months
season_mean <- calc(anoms_nodec[[in_season]], mean, na.rm = T)
# season_mean <- st_as_stars(rotate(season_mean))
return(season_mean)
})
# Plotting the seasons
seas_anom_maps <- function(season, lab_years){
# Center the color scale with Dynamic Limits
temp_limits <- max(abs(values(season_stacks[[season]])), na.rm = T) * c(-1,1)
# Make map
seas_map <- ggplot() +
geom_stars(data = st_as_stars(rotate(season_stacks[[season]]))) +
geom_sf(data = world, fill = "gray30", color = "white", size = .25) +
scale_fill_distiller(palette = "RdBu", na.value = "transparent",
limit = temp_limits, oob = scales::squish) +
guides("fill" = guide_colorbar(title = sst_lab,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(3, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
coord_sf(expand = F) +
map_theme +
labs(title = str_c(season, " - ", lab_years))
return(seas_map)
}
seas_anom_maps("Winter", "2020-2021")
seas_anom_maps("Spring", "2021")
seas_anom_maps("Summer", "2021")
#### Making Maps of Each
seas_anom_maps("Fall", "2021")
# Create an expanded area to see regional patterns beyond the explicit region used
# # Expand the area out to see the larger patterns
# crop_x <- crop_x + c(-2.5, 2.5)
# crop_y <- crop_y + c(-1.5, 1)
# Make a new extent for cropping
region_extent_expanded <- st_sfc(st_polygon(list(
rbind(c(crop_x[[1]], crop_y[[1]]),
c(crop_x[[1]], crop_y[[2]]),
c(crop_x[[2]], crop_y[[2]]),
c(crop_x[[2]], crop_y[[1]]),
c(crop_x[[1]], crop_y[[1]])))))
region_extent_expanded <- st_as_sf(region_extent_expanded)
# Mask the annual average
reg_ann_anom <- mask_nc(ras_obj = ann_anom_ras, mask_shape = region_extent_expanded)
# Center temperature limits
temp_limits <- max(abs(values(reg_ann_anom)), na.rm = T) * c(-1,1) # Dynamic Limits
# Build Map
ggplot() +
geom_stars(data = st_as_stars(reg_ann_anom)) +
geom_sf(data = new_england, fill = "gray30", color = "white", size = .25) +
geom_sf(data = canada, fill = "gray30", color = "white", size = .25) +
geom_sf(data = greenland, fill = "gray30", color = "white", size = .25) +
scale_fill_distiller(palette = "RdBu", na.value = "transparent",
limit = temp_limits, oob = scales::squish) +
guides("fill" = guide_colorbar(title = sst_lab,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(3, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(xlim = crop_x,
ylim = crop_y, expand = F) +
labs(title = str_c("Regional Temperature Anomalies - ", nc_year))
# Mask the Seasons
seasons_masked <- map(season_stacks, mask_nc, region_extent_expanded)
# Plot the seasons
plot_masked_season <- function(season, year_lab){
# Grab Season
reg_seas_anom <- seasons_masked[[season]]
# Get temp limits
temp_limits <- max(abs(values(reg_seas_anom)), na.rm = T) * c(-1,1) # Dynamic Limits
# Build Map
seas_map <- ggplot() +
geom_stars(data = st_as_stars(reg_seas_anom)) +
geom_sf(data = new_england, fill = "gray30", color = "white", size = .25) +
geom_sf(data = canada, fill = "gray30", color = "white", size = .25) +
geom_sf(data = greenland, fill = "gray30", color = "white", size = .25) +
scale_fill_distiller(palette = "RdBu", na.value = "transparent",
limit = temp_limits, oob = scales::squish) +
guides("fill" = guide_colorbar(title = sst_lab,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(3, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(xlim = crop_x,
ylim = crop_y, expand = F) +
labs(title = str_c(season, " - ", year_lab))
return(seas_map)
}
plot_masked_season("Winter", year_lab = "2020-2021")
plot_masked_season("Spring", year_lab = "2021")
plot_masked_season("Summer", year_lab = "2021")
plot_masked_season("Fall", year_lab = "2021")
Looking specifically at the last heatwave event, we can step through how the event progressed over time, and developing pockets or warmer/colder water masses.
# Identify the last heatwave event that happened
last_event <- max(region_heatwaves$mhw_event_no, na.rm = T)
# Pull the dates of the most recent heatwave
last_event_dates <- region_heatwaves %>%
filter(mhw_event_no == last_event) %>%
pull(time)
# Buffer the dates, start 7 days before
event_start <- (min(last_event_dates) - 7)
event_stop <- max(last_event_dates)
date_seq <- seq.Date(from = event_start,
to = event_stop,
by = 1)
# Expand the area out to see the larger patterns
crop_x <- crop_x + c(-2.5, 2.5)
crop_y <- crop_y + c(-1.5, 1)
# Load the heatwave dates
data_window <- data.frame(
time = c(min(date_seq) , max(date_seq) ),
lon = crop_x,
lat = crop_y)
# Pull data off box
hw_stack <- oisst_window_load(data_window = data_window,
anomalies = T, mac_os = "mojave")
#drop any empty years that bug in
hw_stack <- hw_stack[map(hw_stack, class) != "character"]
##### Format the layers and loop through the maps ####
# Grab only current year, format dates
this_yr <- stack(hw_stack)
day_count <- length(names(this_yr))
day_labs <- str_replace_all(names(this_yr), "[.]","-")
day_labs <- str_replace_all(day_labs, "X", "")
day_count <- c(1:day_count) %>% setNames(day_labs)
# Progress through daily timeline to indicate heatwave status and severity
hw_timeline <- region_heatwaves %>%
filter(time %in% as.Date(day_labs))
#### Plot Settings:
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- c(max(values(this_yr), na.rm = T) * -1,
max(values(this_yr), na.rm = T) )
# Plot Heatwave 1 day at a time as a GIF
day_plots <- imap(day_count, function(date_index, date_label) {
# grab dates
heatwaves_st <- st_as_stars(this_yr[[date_index]])
#### 1. Map the Anomalies in Space
day_plot <- ggplot() +
geom_stars(data = heatwaves_st) +
geom_sf(data = new_england, fill = "gray30", color = "white", size = .25) +
geom_sf(data = canada, fill = "gray30", color = "white", size = .25) +
geom_sf(data = greenland, fill = "gray30", color = "white", size = .25) +
geom_sf(data = region_extent,
color = gmri_cols("gmri blue"),
linetype = 2, size = 1,
fill = "transparent") +
scale_fill_distiller(palette = "RdYlBu",
na.value = "transparent",
limit = limit) +
map_theme +
coord_sf(xlim = crop_x,
ylim = crop_y,
expand = T) +
guides("fill" = guide_colorbar(
title = "Sea Surface Temperature Anomaly \u00b0C",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(3, "in"),
frame.colour = "black",
ticks.colour = "black"))
#### 2. Plot the day and the overall anomaly to track dates
date_timeline <- ggplot(data = hw_timeline, aes(x = time)) +
geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
geom_line(aes(y = hwe, color = "Heatwave Event")) +
geom_line(aes(y = cse, color = "Cold Spell Event")) +
geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 3, size = .5) +
geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
scale_color_manual(values = color_vals) +
# Animated Point / line
geom_point(
data = filter(hw_timeline, time == as.Date(date_label)),
aes(time, sst, shape = factor(mhw_event)),
color = gmri_cols("gmri blue"),
size = 3, show.legend = FALSE) +
geom_vline(data = filter(hw_timeline, time == as.Date(date_label)),
aes(xintercept = time),
color = "gray50",
size = 0.5,
linetype = 3,
alpha = 0.8) +
labs(x = "",
y = "",
color = "",
subtitle = "Regional Temperature \u00b0C",
shape = "Heatwave Event") +
theme(legend.position = "bottom")
#### 3. Assemble plot(s)
p_layout <- c(
area(t = 1, l = 1, b = 2, r = 8),
area(t = 3, l = 1, b = 8, r = 8))
# plot_agg <- (date_timeline / day_plot) + plot_layout(heights = c(1, 3))
plot_agg <- date_timeline + day_plot + plot_layout(design = p_layout)
return(plot_agg )
})
walk(day_plots, print)